代码仓库：https://gitee.com/gislxz/zonal-statistics

原文来自个人博客：https://www.gislxz.com/index.php/2023/04/19/zonalstatistics/

个人博客中还有GDAL学习笔记以及GIS相关的学习&开发记录，感兴趣的小伙伴们可以看看：https://www.gislxz.com/

背景

实现这个算法的想法萌生于去年，我的导师指导我去协助参与重庆农业项目，任务是更新重庆市所有地块一年的NDVI数据，这个任务看似简单但工作量极大，除去数据准备外，该任务最核心的步骤就是使用分区统计将栅格格式的NDVI数据赋值到以多边形为组织形式的地块数据中。然而就是这样简单的分区统计功能却耗费了我们大量时间，ArcGIS和QGIS软件提供的分区统计功能并不足以支持我们高效精准的完成这项任务，只能退而求其次选择了中心点赋值法完成了NDVI赋值的任务，但这样不利于数据不确定性的控制，进而影响后续计算与预测，造成误差的传递和放大。

现有GIS软件分区统计功能的不足

Arcgis：因为我只有arcgis10.8的版本，不太清楚新的Arcgis Pro 3的算法实现怎么样，但10.8版本中分区统计算法太慢了，比QGIS慢数十倍，输出是个表格，还要再与原文件进行连接。

QGIS：QGIS处理倒挺快的，问题是不能写入原文件，只能输出一个新矢量文件，那我一次运算12个月的栅格图像，没法批处理运行了。

思来想去还是自己重写一个。

使用GDALWarp来实现

一开始倒是想用现成的API来实现的，GDALWarp中有图像裁切的API可以拿来用。分区统计本质上就是统计每一个多边形内部像元点的值，其实也就是要对图像进行裁切。只不过我们不需要将裁切后的图像保存为新文件，只用把像元值存储的二维数组输入到内存里在对数组进行运算就好了。遗憾的是GDAL学习资料太少，李民录先生写的《GDAL源码剖析与开发指南》介绍了Warp的基本用法，但是并没说怎么输入到内存，只说了输入到内存是用WarpRegionToBuffer这个执行语句。我也是进行了一番尝试，这个Warp注意事项主要有两点，一是输入的AOI范围需要是多边形或多多边形，其坐标应该是像素行列号而不是地理坐标，且要求是WKTGeo格式，好在OGRGeometry类是有exportToWkt的函数。第二点就是使用Warp做裁切还需要设置坐标转换的参数，书中给的示例是用的GDALCreateGenImgProjTransformer2这个函数，输入输入和输出图像对应的GDALDataset就行了，但我们不需要输出到文件，自然也没有创建输出栅格图像的GDALDataset，如果创建势必会极大影响处理速度，于是我去查了其他几个GDALCreateGenImgProjTransformer函数，发现可以这样替代。

调用GDALWarp进行裁剪的完整函数代码

测试后发现调用GDALWarp来实现存在以下几点不足：

1 处理速度比QGIS还是慢四到五倍

2 不能精确捕获NODATA值

3 面积过小的地块发生错误

于是我只能尝试从头实现，自己来写裁切功能。

算法设计思路

在判断多边形与栅格图像相交范围中，我采用了中心点判断法，及像元中心点被多边形所包含，即认为多边形与栅格图像的相交范围包含该像元。

计算多边形与栅格图像相交范围步骤：

1通过栅格图像行列号转地理坐标的转换六参数反推地理坐标转行列号的六参数。（需注意转换六参数中的起始点为左上角像素点的左上角位置而非中心点，因此反算时需要将起始点的坐标分别加上对应方向像素宽度的一半）

2对任一多边形遍历其外环（ExteriorRing），将得到的点坐标转为行列号，通过前后两个点连成的直线计算跨越行号时与对应直线的交点。

3对每一行的交点进行排序，从小到大即从左到右，这样包含在俩点间的像素点必然都在多边形内。对于凹多边形，多多边形，内环等情况该算法都适用。

在获得多边形包含的像元点后即可读取栅格将对应像素点的数值赋给多边形，再进行统计计算。这里我注意到如果遍历多边形以其四至范围去读取栅格并赋值会造成极大的硬盘IO开销，因而将策略调整成逐行读取栅格图像，并对该行所有相交的多边形进行赋值处理。

回忆NDVI赋值的任务，每月对应的一幅NDVI图像的尺寸，角度和范围可以认为是相同的，这就意味着每个多边形计算的相交（掩膜）范围是可以重复利用的，因此我又增加了一个批处理加速的策略，当输入多幅图像时可以使用批处理加速，在处理完每个地块时不释放相交范围和记录像元值缓冲区而只是重置有效像元计数，用更多的内存占用换取处理时间上的减少。

算法结构设计

说明：

1整个代码使用Entrance函数作为入口，打开第一幅栅格图像以获取栅格数据类型从而创建相应ZonalStatistics（模板类），并把打开栅格图像所用的GDALDatasetUniquePtr转给（std::move）ZonalStatistics。

2 ZonalStatistics类中具有一个ParcelManager指针，在构造函数中会构造一个其独有的ParcelManager。ZonalStatistics类包含了栅格图像文件的指针，需要添加的属性信息，其成员函数void poiStatistics()和void MyStatistics()分别实现了借助GDALWrap裁剪图像和底层构建的分区统计算法。

3 ParcelManager作为管理ParcelUnit的类，其关键成员在于具有两个std::forward_list**成员分别用于记录栅格图像某一行需要启动和结算的地块对象，将其加入处理队列或从处理队列中删除并进行结算（批处理加速中只重置计数）。处理队列使用std::unordered_set*结构实现方便插入和删除。

4 ParcelUnit作为地块单元，使用自定义结构体MyExtent记录几何体范围，以便插入ParcelManager的启动和结算队列；使用std::vector**记录相交（掩膜）范围；使用double数组记录包含的像元值。同时其具有一个中心点坐标和对应像元值（如果在栅格图像范围内）以便处理地块过小不包含任一像元中心点的情况。

结果对比

实验数据为

·重庆市2020年10月及2020年11月NDVI合成图像（单波段，尺寸54615×44970，单幅数据大小9.89GB）

·重庆市綦江区耕地地块数据，共761496个几何体，包含凹多边形，多多边形，内环等各种情况。

1本算法与QGIS分区统计算法计算结果基本相同，QGIS对小地块有特殊处理策略，因而造成个别地块结果存在些许差异

2在处理单幅图像中本算法与QGIS处理速度基本一致，在批量处理中，QGIS因为属性表体积的增长而减速，本算法在向属性表中添加列步骤中耗时没有明显变化且统计计算步骤由于优化策略有明显提速。

这次分区功能的实现相当于我学习GDAL过程中第一次的实践尝试，最多的感觉就是参考资料不足，很多时候要去GDAL官方文档里扒人家API的说明。另外对于算法的实现思路也是拍脑袋想的，没有去查相关资料（其实查了一下中文的论文没查到，懒得查英文的）。最近在阅读《计算几何：算法与应用》，才发现原来从上而下的平面扫描是这么经典的思路，还有这么多可以改进的地方，以后还是要多看书多读论文，不能拍脑袋想当然的去写了，否则遇到一些更加复杂困难的功能想出的实现思路一定不是最优解，可能计算复杂度差一俩个数量级。